Accurate interaction energies from perturbation theory based on Kohn-Sham model 



in 
o 
o 

(N 

c 
a 
>—» 



43 

Oh- 
i ' 

S 

1) ■ 
43 ' 

q ; 

o ■ 

<Z3 , 
>v 
43 ' 
Oh: 



> 

(N 

o 

o 
in 
o 

o 

43 

Oh. 



Rafal Podeszwa and Krzysztof Szalewicz 
Department of Physics and Astronomy, University of Delaware, Newark, DE 19716 

(Dated: February 2, 2008) 

The density-functional based symmetry-adapted perturbation theory [SAPT(DFT)] has been ap- 
plied to the argon, krypton, and benzene dimers. It is shown that — at a small fraction of computa- 
tional costs — SAPT(DFT) can provide similar accuracies for the interaction energies as high-level 
wave-function based methods with extrapolations to the complete basis set limits. This accuracy is 
significantly higher than that of any other DFT or DFT-based approaches proposed to date. 
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The intermolecular forces — sometimes called van der 
Waals (vdW) interactions or forces — determine the struc- 
ture and properties of most clusters, liquids, and solids. 
These forces also govern many life processes, such as the 
genetic code replication, protein structure and dynamics, 
and enzymatic actions. Thus, the ability to computation- 
ally predict van der Waals forces is significant for the un- 
derstanding of all these systems. However, although the 
standard wave-function based electronic structure meth- 
ods can in principle be used for such predictions, in prac- 
tice these methods are too time consuming to be ap- 
plied to most systems of interest in biology, even with 
extensive use of the current computers capabilities. The 
density-functional theory (DFT) methods are much less 
time-consuming, however, the currently existing versions 
of DFT fail to describe the dispersion interaction, an im- 
portant part of the van der Waals force. This problem is 
due to the fact that the dispersion forces result from long- 
range correlations between electrons, whereas the current 
exchange-correlation potentials model only short-range 
correlation effects. Many authors add the asymptotic ex- 
pansion of the dispersion energy to the DFT interaction 
energies, which inherently includes some double counting 
(see Ref.Q for a discussion of these issues). 

Occasionally, for a specific system, one of the variants 
of DFT can give reasonably good predictions of inter- 
action energies, at least in some regions of a potential 
energy surface. This fact encouraged some authors to 
build system-specific potentials fitted to a number of grid 
points on the potential energy surface computed using 
a wave-function based method. For example, Boese et 
al. optimized an ammonia-specific potential by ad- 
justing some parameters in one of the standard function- 
als. Recently, Lilienfeld et al. 0] proposed to use atom- 
centered nonlocal effective core potentials with parame- 
ters adjusted for specific systems. These methods do not 
offer physically motivated improvement of the density- 
functional formalism but rather rely on cancellations of 
errors. Also, a number of wave-function calculations are 
needed to optimize the parameters in the functionals, 
which limits the range of applications to systems that 
can be treated with the latter methods (unless the pa- 
rameters can be shown to be transferable). 



An approach which uses the specific characteristic of 
the dispersion interaction has recently been presented by 
Dion et al. Q . The method was denoted by the authors 
as vdW-DF. It adds a nonlocal correlation energy part 
to existing functionals. This term models the dispersion 
energy utilizing approximate density response functions. 
The method predicts interaction energies of the systems 
investigated in Ref. qualitatively (to within a factor of 
about 1.5, see below). 

Another approach to the calculations of interaction 
energies for large molecules was developed by Mis- 
quitta et al. la, ||| and independently by Hesselmann 
and Jansen Q, following ideas of Williams and Cha- 
balowski ||. This approach is based on symmetry- 
adapted perturbation theory (SAPT) Q, but utilizes 
the description of the interacting monomers in terms 
of Kohn-Sham (KS) orbitals, orbital energies, and 
frequency-dependent density susceptibility (FDDS) func- 
tions. The DFT-based SAPT will be called SAPT(DFT). 
This method can be shown to be potentially exact for all 
major components of the interaction energy (asymptot- 
ically for exchange interactions) in the sense that these 
components would be exact if the DFT description of the 
monomers were exact [fulfil I lflj . Applications to a number 
of small dimers have shown that SAPT(DFT) provides 
surprisingly accurate results, sometimes more accurate 
than the standard SAPT at the currently programmed 
level [H Ej. 

The regular SAPT method involves expansions in pow- 
ers of the intermonomer interaction operator V and 
the intramonomer correlation operator W, the so-called 
M0ller-Plesset (MP) potential. The terms proportional 
to powers of W describe the effects of intramonomer 
electron correlation on the interaction energy. Simi- 
larly as in the electronic-structure many-body pertur- 
bation theory (MBPT) or coupled-cluster (CC) meth- 
ods, these terms are expensive to compute, with CPU 
times scaling as a high power of system size N — the sev- 
enth power if the complete currently programmed set 
of SAPT components is included. This scaling is the 
same as for MBPT in the fourth order (MP4) or CC 
including single, double, and noniterated triple excita- 
tions [CCSD(T)]. In SAPT(DFT), no terms of this type 
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appear, as the intramonomer correlation effects are ac- 
counted for by DFT. Therefore, SAPT(DFT) scales as 
only A^ 5 , i.e., a SAPT(DFT) calculation is generally or- 
ders of magnitude faster than a regular SAPT calculation 
at the complete currently programmed level. This com- 
putational advantage is further significantly increased by 
the superior basis set convergence of SAPT(DFT) com- 
pared to the wave-function based electron correlation 
methods. The latter methods converge slowly due to 
the necessity of reproducing the intramonomer electron- 
electron cusps by expansions in products of one-electron 
functions. Such expansions do not appear in DFT. In 
SAPT (DFT), the orbital-product expansions are present 
in the expressions for the dispersion energy, however, it 
has been shown [1(1 IllT ] that this component (similarly 
as the SAPT dispersion energy of zeroth order in W) 
can be saturated in reasonably small basis sets provided 
that "midbond" functions are used, i.e., basis functions 
are placed at a point between the two monomers. It 
appears that, in many cases, polarized triple-zeta (TZ) 
quality bases give SAPT(DFT) interaction energy com- 
ponents converged to a similar number of digits as the 
regular SAPT components in polarized quadruple-zeta 
(QZ) bases. This results in a difference in the basis set 
size of about a factor of two, a ratio of 2 4 in computer 
time at the MP4 or CCSD(T) level. Our current imple- 
mentation of SAPT(DFT) is using an interface to the 
time-dependent DFT (TD-DFT) part CADPAC G3 to 
compute FDDS's, which is the time- limiting step of the 
calculations. An optimized TD-DFT program now un- 
der development should decrease time requirements of 
this step by at least one order of magnitude. This will 
make the largest calculations described here comparable 
to the supermolecular DFT calculations in the same ba- 
sis. Thus, the SAPT(DFT) method is not prohibitively 
expensive, as stated by the authors of Ref. 0, and it 
has already been applied t o sy stems containing about 
40 atoms and 200 electrons [l3| . 

We present here the first application of SAPT(DFT) 
to relatively large interacting monomers. In order to 
compare with the vdW-DF method, we have chosen the 
same systems as investigated in Ref. 0: Ar2, Kr 2 , and 
the benzene dimer. We have used the following Carte- 
sian basis sets: aug-cc-pV5Z Q for Ar 2 (with g and 
h functions removed due to restrictions of CADPAC), 
aug-cc-pVTZ for Kr 2 , and a polarized double-zeta 
size basis with polarization coefficients optimized on dis- 
persion energy ^| for benzene. In all cases, we used 
a set of midbond functions consisting of s and p func- 
tions with exponents 0.9, 0.3, 0.1, and d and / functions 
with exponents 0.6 and 0.2, placed at the center of mass 
of each dimer. We used the monomer-centered "plus" 
form of basis sets 0] for the argon (including d functions 
on the interacting partner) and benzene dimers, and the 
full dimer-centered form Kr 2 . For the benzene dimer, we 
considered the parallel "sandwich" configuration, and the 
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FIG. 1: (Color online). Interaction energy of the ar gon 
dimer. SAPT(DFT) — this work. CCSD(T]/CBS— Ref. M 
Benchmark— Aziz, Ref-HJ vdW-DF— Ref . |1 



monomer geometry was taken from Ref. We employed 
PBE0 ^lSg DFT functional with the asymptotic correc- 
tion [f| for all the systems. In addition, the B97-2 0] 
functional was applied at near minimum geometries for 
all dimers and at all points for the argon dimer. Effects 
of the third and higher orders in V have been neglected. 
In TD-DFT, we have used the standard PBE0 or B97-2 
kernels for argon and krypton, and the LDA kernel for 
benzene. The use of the LDA kernel provides a consid- 
erable speedup of calculations. To check the accuracy of 
this approximation, we performed a single-point calcula- 
tion using the PBE0 kernel for the benzene dimer. The 
error in the dispersion energy was smaller than 1%. 

Figure ^ presents the interaction potential for the ar- 
gon dimer. The benchmark results are the empirical 
potential of Aziz j^] and the CCSD(T) potential with 
extrapolation to the complete basis set (CBS) limit by 
Patkowski et al. [2(J. The two curves are almost indis- 
tinguishable, showing the very high level of agreement 
between ab initio theory and experiment for this sys- 
tem. The SAPT (DFT) calculations are very close to the 
benchmarks, within about 2 cm -1 or 2% at the mini- 
mum geometry. In contrast, the vdW-DF method gives 
a curve which is about 1.6 times too deep, and the min- 
imum position is shifted by 0.2 A. 

Our results for the krypton dimer are displayed in 
Fig. [21 For this system, the benchmark curve is given by 
the empirical fit of Dham et al. H^. The CCSD(T) curve 
computed by Slavicek et al. [2^ represents the best litera- 
ture non-relativistic theoretical result and is in fact at the 
limit of what the ab initio methods are capable of achiev- 
ing at the present time. The SAPT(DFT) curve agrees 
with the benchmark slightly better than the CCSD(T) 
one, and much better than the curve produced by the 
vdW-DF method The latter curve is about a factor 
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FIG. 2: (Color online). Interaction energy of the kryp- 
ton dimer. SAPT(DFT)— this work. Benchmark— Dham et 
al. |22^ . CCSD(T) — calculations in aug-cc-pV5Z-|-,spd/'</ basis 
by Slavicek et al. vdW-DF— Ref.f" 
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FIG. 3: (Color online). Interaction energy between two 
benzene molecules in the parallel "sandwich" configuration. 
SAPT(DFT) — this work. CCSD(T)— Ref. E3 (Model II). 
MP2 Tsuzuki et al. computed in the aug(d,p)-6-311G** 

basis. MP2 Sinnokrot et al . [24l| — computed in the aug-cc- 
pVTZ basis. vdW-DF— Ref.liDFT Lilienfeld et al— Ref.0 



of 1.4 too deep and has the minimum position shifted by 
about 0.2 A. 

Finally, the benzene dimer results are presented in 
Fig. El For this system, there is no highly reliable bench- 
mark available. The best ab initio calculations for a 
number of monomer separations are those by Tsuzuki 
et al. 17j at the CCSD(T) level. The values displayed 
were termed as "Model II" and included the MP2 ener- 
gies computed in a polarized DZ-quality basis augmented 
with a single diffuse polarization set [aug(d,p)-6-311G**] 
and the ACCSD(T) = CCSD(T)-MP2 component com- 
puted in a DZ basis polarized only on carbons. Below 



TABLE I: Dimer interaction energies near minima. 
SAPT(DFT) results obtained using the PBE0 and B97-2 
functionals are compared with benchmark results. For the 
argon and krypton dimers the energies are in cm -1 , for the 
benzene dimer in kcal/mol. 





R(k) 


PBE0 


B97-2 


Benchmark 


Ar 2 


3.75 


-97.46 


-97.81 


-99.64" 


Kr 2 


4.00 


-140.49 


-141.36 


-139.87 6 


(C6H6)2 


3.80 


-1.48 


-1.52 


-1.48 c , -1.81 d 



"Ref. 
6 Ref. 
c Ref. 
d Ref. 



-CCSD(T) (Model III). 
-CCSD(T). 



we will quote the MP2 and ACCSD(T) values in paren- 
theses. Model II gives -1.62 (-2.85, 1.24) kcal/mol at 
3.8 A. The authors computed also an extrapolated inter- 
action energy ("Model III", only at 3.8 A) by perform- 
ing a CBS extrapolation at the MP2 level and scaling 
the ACCSD(T) component, which gave -1.48 (-3.28, 
1.80) kcal/mol. Accurate calculations for the benzene 
dimer were also performed by Sinnokrot et al. |24j |. but 
only at a single point, R = 3.7 A. These authors com- 
puted the MP2 energies in bases up to aug-cc-pVQZ and 
ACCSD(T) in aug-cc-pVDZ, which gives -1.54 (-3.37, 
1.83) kcal/mol. If one extrapolates the MP2 TZ-QZ re- 
sults of Ref. using the X~ 3 extrapolation scheme, one 
obtains the value of —3.45 kcal/mol, not much different 
from the best calculated result. Sinnokrot et al. com- 
puted the MP2 energy using also the so-called MP2-R12 
explicitly correlated basis and their recommended inter- 
action energy is —1.81 (—3.64, 1.83) kcal/mol. The fairly 
large discrepancy between the extrapolated and MP2- 
R12 results could be due to an insufficient convergence 
of the resolution of identity applied in the MP2-R12 ap- 
proach in the [spdf/spd] bases. Based on these consid- 
erations, one has to assume that the uncertainty of the 
CCSD(T) curve in Figgis about ±0.2 kcal/mol. 

Figure shows that the agreement of SAPT(DFT) 
with the CCSD(T) results is very good, in particular tak- 
ing into account the uncertainty of the latter. In fact, in 
view of this uncertainty, the SAPT(DFT) results provide 
an independent set of the best current estimates of the 
exact interaction energies for the benzene dimer. The 
vdW-DF curve is deeper by about a factor of 1.4 and the 
minimum is shifted by about 0.2 A. Figure|3includes also 
the MP 2 results from Ref. El the DFT results of Lilien- 
feld et al. from Ref. 0, and the MP2 results of Sinnokrot 
et al. 2^ which were used to calibrate the Lilienfeld et 
al. DFT functional. The differences between the two 
MP 2 curves are consistent with the sizes of the basis set 
effects discussed above. For the benzene dimer, the MP2 
level of theory is not adequate and therefore the results of 
Lilienfeld et al. are very far from the CCSD(T) bench- 
mark. If these authors had chosen to fit their functional 
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TABLE II: Interaction energy components (in kcal/mol) for 
the benzene dimer at the parallel "sandwich" geometry ob- 
tained with SAPT(DFT) compared to those obtained with the 
regular SAPT and a modified aug-cc-pVDZ basis in Ref. l2rj 

PBEO/3.8 A PBEO/3.7 A SAPT/3.7 A" 



electrostatic 


0.09 


-0.28 


-0.97 


1st order exchange 


3.66 


4.94 


6.03 


induction 


-1.24 


-1.69 


-2.14 


exchange-induction 


1.02 


1.46 


1.95 


dispersion 


-5.51 


-6.44 


-7.47 


exchange-dispersion 


0.50 


0.65 


0.94 


total 


-1.48 


-1.37 


-1.66 



"See Ref. for the SAPT components included at each level. 



to the CCSD(T) level of theory, the calculations would 
have become significantly more expensive. This empha- 
sizes the fact that the DFT results of Ref. can only be 
as accurate as the underlying wave-function based calcu- 
lations. Even within the MP2 model, Fig. [3] shows that 
the agreement between the DFT and MP2 energies dete- 
riorates quickly in the regions farther from the minimum. 

It should be emphasized that, in contrast to the su- 
permolecular DFT approach, SAPT(DFT) is relatively 
insensitive to the choice of the DFT functional. Sev- 
eral published papers show that the results given by the 
former method with different choices of the functionals 
can differ by an order of magnitude and can be of the 
wrong sign (see, e.g., Refs. OHHj). The results displayed 
in Fig. [Hand in Table [Q show that the SAPT(DFT) re- 
sults computed with PBE0 and B97-2, two functionals 
developed using very different principles, agree very well, 
with discrepancies being of the same magnitude as other 
uncertainties of the SAPT(DFT) approach. 

An advantage of the SAPT approach is that it gives di- 
rectly the physical components of the interaction energy. 
In Tabled we present such components for the benzene 
dimer and compare with SAPT calculations by Sinnokrot 
and Sherrill [26|. We have omitted in this comparison 
terms beyond the second order in V computed in RefErj 
These terms amount to —0.14 kcal/mol. Generally, the 
agreement between SAPT (DFT) and SAPT is good tak- 
ing into account different levels of intramonomer corre- 
lation effects, different basis sets, and different monomer 
geometries used. 

In summary, we have shown that SAPT(DFT) is capa- 
ble to achieve the accuracy of intermolecular interaction 
energies similar to that of the CCSD(T)/CBS approach 
at a very small fraction of computational costs. This ac- 
curacy is much higher than that of the vdW-DF method 
of Dion et al. 4]. The SAPT(DFT) method is rigor- 
ously valid for all separations between the interacting 
molecules. It therefore constitutes the most accurate cur- 
rent approach for practical calculations of interactions for 
large monomers, containing as many as 20 atoms each. 
This development will bring some important biophysi- 



cal applications within reach of computational physics 
(e.g., interactions involving DNA bases, small polypep- 
tides, and sugars). 

The authors are grateful to Bogumil Jeziorski for read- 
ing the manuscript and for valuable advice. This research 
was supported by a grant from ARO. 
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